function H_tp1_peers = mean_peers(H1,H2,PS1,PS2, gamma,fp,set_kids2,fp_parfor,t) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Simulate Realization of Peer Quality;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

%Child Utility Function (equation 15);
U_links1 = gamma(1).*ones(size(H1))  +  ...
           gamma(2).*log(H1) + gamma(3).*log(H2) + gamma(4).*(log(H1) - log(H2)).^(2) + ...
           gamma(5).*( (log(H1) - log(H2)).^(2) ).*(log(H1) - log(H2)>0).*PS1;                                                 
       
U_links2 = gamma(1).*ones(size(H1))  +  ...
           gamma(2).*log(H2) + gamma(3).*log(H1) + gamma(4).*(log(H1) - log(H2)).^(2) + ...
           gamma(5).*( (log(H1) - log(H2)).^(2) ).*(log(H2) - log(H1)>=0).*PS2;         

%Probability of a friendship (equation 16);
Prob_Links_temp = (1./(1+1./exp(U_links1))).*(1./(1+1./exp(U_links2))) ;    
Prob_Links = repmat( Prob_Links_temp, [1 , 1 , fp.tot_draws_mcintegration_network]) ; 

%Simulate friendships;
Links = ( Prob_Links >= fp_parfor.peers_backward(:,:,:,t));        
%Utility of Friendships;
Child_util = repmat( U_links1, [1 , 1 , fp.tot_draws_mcintegration_network]) ; 
Child_util = Child_util.*Links; %Equation 7;
peers = repmat(set_kids2',[size(Prob_Links,1) ,1, fp.tot_draws_mcintegration_network]) ;         
temp = Links.*peers;
%Collect the various simulations of average peers;
H_prime = nansum(temp,2)./sum(temp~=0 &  isnan(temp)==0,2) ;  
H_prime = reshape(H_prime,[size(H_prime,1),size(H_prime,3)] ) ;                        
        
H_tp1_peers.H_prime = H_prime;
H_tp1_peers.U_links1   = Child_util ;